#if !defined(MTX_FIXED_H)
#define MTX_FIXED_H

#include "mtx/types.h"
#include <cassert>

namespace mtx {

template <size_type n> class fixed_vector_base {
public:
	typedef double * iterator;
	typedef const double * const_iterator;
public:
	enum { _size = n };
	double _data[n];
public:
	fixed_vector_base() 
		{ }
	size_type size() const
		{ return n; }
	const double & operator () (size_type i) const
		{ assert(i < n); return _data[i]; }
	double & operator () (size_type i)
		{ assert(i < n); return _data[i]; }
	const_iterator begin() const
		{ return _data; }
	const_iterator end() const
		{ return _data + n; }
	iterator begin()
		{ return _data; }
	iterator end()
		{ return _data + n; }
};

template <size_type m, size_type n> class fixed_matrix_base {
public:
	typedef double * iterator;
	typedef const double * const_iterator;
public:
	enum { _rows = m, _cols = n, _size = m * n };
	double _data[_cols][_rows];
public:
	fixed_matrix_base() 
		{ }
	size_type rows() const
		{ return _rows; }
	size_type cols() const
		{ return _cols; }
	const double & operator () (size_type i, size_type j) const
		{ assert(i < _rows && j < _cols); return _data[j][i]; }
	double & operator () (size_type i, size_type j)
		{ assert(i < _rows && j < _cols); return _data[j][i]; }
	const_iterator begin() const
		{ return reinterpret_cast<const_iterator>(_data); }
	const_iterator end() const
		{ return reinterpret_cast<const_iterator>(_data) + _size; }
	iterator begin()
		{ return reinterpret_cast<iterator>(_data); }
	iterator end()
		{ return reinterpret_cast<iterator>(_data) + _size; }
};

template <size_type m, size_type n = 1>
struct fixed {
public:
	typedef vector< fixed_vector_base<m> > vector;
	typedef matrix< fixed_matrix_base<m, n> > matrix;
};			  

template <size_type n>
vector< fixed_vector_base<n> > operator + (const vector< fixed_vector_base<n> > & x, const vector< fixed_vector_base<n> > & y)
{
	vector< fixed_vector_base<n> > t(x);
	t += y;
	return t;
}

template <size_type n>
vector< fixed_vector_base<n> > operator - (const vector< fixed_vector_base<n> > & x, const vector< fixed_vector_base<n> > & y)
{
	vector< fixed_vector_base<n> > t(x);
	t -= y;
	return t;
}

template <size_type n>
vector< fixed_vector_base<n> > operator * (const vector< fixed_vector_base<n> > & v, double x)
{
	vector< fixed_vector_base<n> > t(v);
	t *= x;
	return t;
}

template <size_type n>
vector< fixed_vector_base<n> > operator / (const vector< fixed_vector_base<n> > & v, double x)
{
	vector< fixed_vector_base<n> > t(v);
	t /= x;
	return t;
}

template <size_type m, size_type n>
matrix< fixed_matrix_base<m, n> > operator + (const matrix< fixed_matrix_base<m, n> > & x, const matrix< fixed_matrix_base<m, n> > & y)
{
	matrix< fixed_matrix_base<m, n> > t(x);
	t += y;
	return t;
}

template <size_type m, size_type n>
matrix< fixed_matrix_base<m, n> > operator - (const matrix< fixed_matrix_base<m, n> > & x, const matrix< fixed_matrix_base<m, n> > & y)
{
	matrix< fixed_matrix_base<m, n> > t(x);
	t -= y;
	return t;
}

template <size_type m, size_type n>
matrix< fixed_matrix_base<m, n> > operator * (const matrix< fixed_matrix_base<m, n> > & mat, double x)
{
	matrix< fixed_matrix_base<m, n> > t(mat);
	t *= x;
	return t;
}

template <size_type n>
vector< fixed_vector_base<n> > operator * (const matrix< fixed_matrix_base<n, n> > & m, const vector< fixed_vector_base<n> > & v)
{
	vector< fixed_vector_base<n> > t;
	for (size_type i = 0; i < n; ++i)
		t(i) = m.row(i) * v;
	return t;
}

template <size_type m, size_type n>
matrix< fixed_matrix_base<m, n> > operator / (const matrix< fixed_matrix_base<m, n> > & mat, double x)
{
	matrix< fixed_matrix_base<m, n> > t(mat);
	t /= x;
	return t;
}

template <size_type n>
matrix< fixed_matrix_base<n, n> > operator * (const matrix< fixed_matrix_base<n, n> > & m1, const matrix< fixed_matrix_base<n, n> > & m2)
{
	matrix< fixed_matrix_base<n, n> > t;
	for (size_type i = 0; i < n; ++i)
		for (size_type j = 0; j < n; ++j)
			t(i, j) = m1.row(i) * m2.col(j);
	return t;
}

inline void inverse(const matrix< fixed_matrix_base<3, 3> > & m, matrix< fixed_matrix_base<3, 3> > & _m)
{
	double det = 
		+ m(0, 0) * (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2)) 
		- m(0, 1) * (m(1, 0) * m(2, 2) - m(2, 0) * m(1, 2))
		+ m(0, 2) * (m(1, 0) * m(2, 1) - m(2, 0) * m(1, 1));
	_m(0, 0) = + (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2));
	_m(0, 1) = - (m(1, 0) * m(2, 2) - m(2, 0) * m(1, 2));
	_m(0, 2) = + (m(1, 0) * m(2, 1) - m(2, 0) * m(1, 1));
	_m(1, 0) = - (m(0, 1) * m(2, 2) - m(2, 1) * m(0, 2));
	_m(1, 1) = + (m(0, 0) * m(2, 2) - m(2, 0) * m(0, 2));
	_m(1, 2) = - (m(0, 0) * m(2, 1) - m(2, 0) * m(0, 1));
	_m(2, 0) = + (m(0, 1) * m(1, 2) - m(1, 1) * m(0, 2));
	_m(2, 1) = - (m(0, 0) * m(1, 2) - m(1, 0) * m(0, 2));
	_m(2, 2) = + (m(0, 0) * m(1, 1) - m(1, 0) * m(0, 1));

}

} // mtx


#endif // MTX_FIXED_H
